Load the necessary libraries
library(car) # for regression diagnostics
library(broom) # for tidy output
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(brms)
library(tidybayes)
library(bayesplot)
library(broom.mixed)
library(rstan)
library(patchwork)
library(DHARMa)
library(standist) # for visualizing distributions
library(rstanarm)
library(ggeffects)
library(DHARMa)
Some ornithologists were interested in the degree of sibling negotiations in owl chicks. Specifically, they wanted to explore how sibling negotiations were affected by feeding satiety and the sex of the parent returning to the nest. The ornithologists had accessed to a number of owl nests and were able to count (via recording equipment) the number of sibling negotiations (calls) that the owl chicks made when the parent returned to the nest.
We could hypothesise that the chicks might call more if they were hungry. As part of the investigation, the researchers were able to provided supplementary food. As such, they were able to manipulate the conditions such that sometimes the chicks in a nest would be considered deprived of supplementary food and at other times they were satiated.
As a parent returned, the researchers recorded the number of sibling negotiations (calls) along with the sex of the parent. Since the number of calls is likely to be a function of the number of chicks (the more chicks the more calls), the researchers also counted the number of siblings in the brood.
Each nest was measured on multiple occasions. Hence, we must include the nest as a random effect to account for the lack of independence between observations on the same set of siblings.
owls <- read_csv("../public/data/owls.csv", trim_ws = TRUE)
owls %>% glimpse()
## Rows: 599
## Columns: 7
## $ Nest <chr> "AutavauxTV", "AutavauxTV", "AutavauxTV", "Autavaux…
## $ FoodTreatment <chr> "Deprived", "Satiated", "Deprived", "Deprived", "De…
## $ SexParent <chr> "Male", "Male", "Male", "Male", "Male", "Male", "Ma…
## $ ArrivalTime <dbl> 22.25, 22.38, 22.53, 22.56, 22.61, 22.65, 22.76, 22…
## $ SiblingNegotiation <dbl> 4, 0, 2, 2, 2, 2, 18, 4, 18, 0, 0, 3, 0, 3, 3, 6, 0…
## $ BroodSize <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, …
## $ NegPerChick <dbl> 0.8, 0.0, 0.4, 0.4, 0.4, 0.4, 3.6, 0.8, 3.6, 0.0, 0…
Let start by declaring the categorical variables and random effect as factors.
## Amount of Sibling negotiation (vocalizations when parents are absent)
## Foot treatment (deprived or satiated
## Sex of parent
## Arrival time of parent
## Nest as random
## Brood size offset
owls <- owls %>% mutate(
Nest = factor(Nest),
FoodTreatment = factor(FoodTreatment),
SexParent = factor(SexParent),
NCalls = SiblingNegotiation
)
As the response represents counts (the number of calls), it would make sense to start by considering a Poisson model. We could attempt to model the response as the number of calls divided by the brood size, but this would result in a response that has no natural distribution.
Instead, if we include brood size as an offset, it will standardise the effects according to brood size (similar to having divided the response by brood size), yet retain the Poisson nature of the response.
The effects of offsets, unlike regular covariates, are not estimated. Rather they are assumed to be 1 (on the link scale). This means that since Poisson uses a log link, then the offset should be of a logged version of the brood size.
Model formula: \[ y_i \sim{} \mathcal{Pois}(\lambda_i)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of food treatment, sex of parent, arrival time (and various interactions) on the number of sibling negotiations. Brood size was also incorporated as an offset. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual nests.
Perhaps we could start off by exploring the main fixed effects. To mimic the log-link, we will use a log-transformed y axis. Since there may well be zeros (no calls detected), we will use a pseudo log scale). We will also include the raw data (jittered and dodged)
ggplot(data = owls, aes(y = NCalls, x = FoodTreatment, color = SexParent)) +
geom_violin() +
geom_point()
ggplot(data = owls, aes(y = NCalls, x = FoodTreatment, color = SexParent)) +
geom_violin() +
geom_point(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.9))
ggplot(data = owls, aes(y = NCalls, x = FoodTreatment, color = SexParent)) +
geom_violin() +
geom_point(position = position_jitterdodge(jitter.height = 0, dodge.width = 1)) +
scale_y_continuous(trans = scales::pseudo_log_trans())
Now, a similar plot separated for each nest.
ggplot(data = owls) +
geom_point(aes(y = NCalls, x = FoodTreatment, color = SexParent), position = position_dodge(0.5)) +
facet_wrap(~Nest)
It might also be worth establishing that there is a linear relationship between the number of calls and brood size. Again, we will mimic the use of the log-link by transforming the axes.
ggplot(data = owls, aes(y = NCalls, x = BroodSize, color = SexParent)) +
geom_point() +
geom_smooth(method = "lm") +
facet_grid(~FoodTreatment) +
scale_y_continuous(trans = scales::pseudo_log_trans()) +
scale_x_log10()
In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over-fitting) and help stabilise the computations.
owls.rstanP <- stan_glmer(NCalls ~ FoodTreatment * SexParent +
offset(log(BroodSize)) + (1 | Nest),
data = owls,
family = poisson(link = "log"),
refresh = 0,
iter = 5000,
warmup = 2000,
thin = 10,
chains = 3,
cores = 3
)
owls.rstanP %>% prior_summary()
## Priors for model '.'
## ------
## Intercept (after predictors centered)
## ~ normal(location = 0, scale = 2.5)
##
## Coefficients
## Specified prior:
## ~ normal(location = [0,0,0], scale = [2.5,2.5,2.5])
## Adjusted prior:
## ~ normal(location = [0,0,0], scale = [5.01,5.08,5.68])
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details
This tells us:
2.5 * sd(owls$NCalls)
## [1] 16.67704
model.matrix(~ FoodTreatment * SexParent, data = owls) %>%
apply(2, sd) %>%
(function(x) 2.5 / x)
## (Intercept) FoodTreatmentSatiated
## Inf 5.007569
## SexParentMale FoodTreatmentSatiated:SexParentMale
## 5.080651 5.679931
Lets now run with priors only so that we can explore the range of values they allow in the posteriors.
owls.rstanarmP1 <- update(owls.rstanP, prior_PD = TRUE)
owls.rstanarmP1 %>%
ggpredict(~ FoodTreatment * SexParent) %>%
plot(add.data = TRUE, jitter = c(0.25, 0)) +
scale_y_continuous("", trans = scales::pseudo_log_trans())
Conclusions:
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):
I will also overlay the raw data for comparison.
owls.rstanarmP2 <- stan_glmer(NCalls ~ FoodTreatment * SexParent +
offset(log(BroodSize)) + (1 | Nest),
data = owls,
family = poisson(link = "log"),
prior_intercept = normal(0, 2.5, autoscale = FALSE),
prior = normal(0, 5, autoscale = FALSE),
prior_covariance = decov(1, 1, 1, 1),
refresh = 0,
iter = 5000,
prior_PD = TRUE,
warmup = 2000,
thin = 10,
chains = 3,
cores = 3
)
owls.rstanarmP2 %>%
ggpredict(~ FoodTreatment * SexParent) %>%
plot(add.data = TRUE, jitter = c(0.25, 0)) +
scale_y_continuous("", trans = scales::pseudo_log_trans())
In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over fitting) and help stabilise the computations.
Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.
owls.form <- bf(NCalls ~ FoodTreatment*SexParent +
offset(log(BroodSize)) + (1|Nest),
family=poisson(link='log'))
options(width=150)
owls.form %>% get_prior(data = owls)
## prior class coef group resp dpar nlpar bound source
## (flat) b default
## (flat) b FoodTreatmentSatiated (vectorized)
## (flat) b FoodTreatmentSatiated:SexParentMale (vectorized)
## (flat) b SexParentMale (vectorized)
## student_t(3, 1.6, 2.5) Intercept default
## student_t(3, 0, 2.5) sd default
## student_t(3, 0, 2.5) sd Nest (vectorized)
## student_t(3, 0, 2.5) sd Intercept Nest (vectorized)
options(width=80)
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):
Note, for hierarchical models, the model will tend to want to have a large \(sigma\) in order to fit the data better. It is a good idea to regularise this tendency by applying a prior that has most mass around zero. Suitable candidates include:
I will also overlay the raw data for comparison.
owls %>%
group_by(FoodTreatment, SexParent) %>%
summarise(
log(median(NCalls)),
log(mad(NCalls))
)
standist::visualize("normal(1.8,5)", xlim = c(0, 20))
standist::visualize("student_t(3, 0, 2.5)",
"cauchy(0,1)",
xlim = c(-10, 25)
)
priors <- prior(normal(1.8, 5), class = "Intercept") +
prior(normal(0, 2), class = "b") +
prior(cauchy(0, 1), class = "sd")
owls.form <- bf(NCalls ~ FoodTreatment * SexParent +
offset(log(BroodSize)) + (1 | Nest),
family = poisson(link = "log")
)
owls.brm2 <- brm(owls.form,
data = owls,
prior = priors,
sample_prior = "yes",
iter = 5000,
warmup = 2500,
chains = 3,
cores = 3,
thin = 10,
refresh = 0
)
owls.form <- bf(NCalls ~ FoodTreatment * SexParent +
offset(log(BroodSize)) + (FoodTreatment * SexParent | Nest),
family = poisson(link = "log")
)
owls.brm3 <- brm(owls.form,
data = owls,
prior = priors,
sample_prior = "yes",
iter = 5000,
warmup = 2500,
chains = 3,
cores = 3,
thin = 10,
refresh = 0,
control = list(adapt_delta = 0.99)
)
(l.1 <- owls.brm2 %>% loo())
##
## Computed from 750 by 599 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2639.3 86.0
## p_loo 148.7 11.7
## looic 5278.7 172.0
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 584 97.5% 119
## (0.5, 0.7] (ok) 11 1.8% 60
## (0.7, 1] (bad) 3 0.5% 11
## (1, Inf) (very bad) 1 0.2% 6
## See help('pareto-k-diagnostic') for details.
(l.2 <- owls.brm3 %>% loo())
##
## Computed from 750 by 599 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2430.8 75.5
## p_loo 301.3 20.8
## looic 4861.6 150.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 507 84.6% 81
## (0.5, 0.7] (ok) 49 8.2% 51
## (0.7, 1] (bad) 37 6.2% 19
## (1, Inf) (very bad) 6 1.0% 1
## See help('pareto-k-diagnostic') for details.
loo_compare(l.1, l.2)
## elpd_diff se_diff
## . 0.0 0.0
## . -208.6 47.6
owls.brm3 %>% get_variables()
## [1] "b_Intercept"
## [2] "b_FoodTreatmentSatiated"
## [3] "b_SexParentMale"
## [4] "b_FoodTreatmentSatiated:SexParentMale"
## [5] "sd_Nest__Intercept"
## [6] "sd_Nest__FoodTreatmentSatiated"
## [7] "sd_Nest__SexParentMale"
## [8] "sd_Nest__FoodTreatmentSatiated:SexParentMale"
## [9] "cor_Nest__Intercept__FoodTreatmentSatiated"
## [10] "cor_Nest__Intercept__SexParentMale"
## [11] "cor_Nest__FoodTreatmentSatiated__SexParentMale"
## [12] "cor_Nest__Intercept__FoodTreatmentSatiated:SexParentMale"
## [13] "cor_Nest__FoodTreatmentSatiated__FoodTreatmentSatiated:SexParentMale"
## [14] "cor_Nest__SexParentMale__FoodTreatmentSatiated:SexParentMale"
## [15] "r_Nest[AutavauxTV,Intercept]"
## [16] "r_Nest[Bochet,Intercept]"
## [17] "r_Nest[Champmartin,Intercept]"
## [18] "r_Nest[ChEsard,Intercept]"
## [19] "r_Nest[Chevroux,Intercept]"
## [20] "r_Nest[CorcellesFavres,Intercept]"
## [21] "r_Nest[Etrabloz,Intercept]"
## [22] "r_Nest[Forel,Intercept]"
## [23] "r_Nest[Franex,Intercept]"
## [24] "r_Nest[GDLV,Intercept]"
## [25] "r_Nest[Gletterens,Intercept]"
## [26] "r_Nest[Henniez,Intercept]"
## [27] "r_Nest[Jeuss,Intercept]"
## [28] "r_Nest[LesPlanches,Intercept]"
## [29] "r_Nest[Lucens,Intercept]"
## [30] "r_Nest[Lully,Intercept]"
## [31] "r_Nest[Marnand,Intercept]"
## [32] "r_Nest[Montet,Intercept]"
## [33] "r_Nest[Murist,Intercept]"
## [34] "r_Nest[Oleyes,Intercept]"
## [35] "r_Nest[Payerne,Intercept]"
## [36] "r_Nest[Rueyes,Intercept]"
## [37] "r_Nest[Seiry,Intercept]"
## [38] "r_Nest[Sevaz,Intercept]"
## [39] "r_Nest[StAubin,Intercept]"
## [40] "r_Nest[Trey,Intercept]"
## [41] "r_Nest[Yvonnand,Intercept]"
## [42] "r_Nest[AutavauxTV,FoodTreatmentSatiated]"
## [43] "r_Nest[Bochet,FoodTreatmentSatiated]"
## [44] "r_Nest[Champmartin,FoodTreatmentSatiated]"
## [45] "r_Nest[ChEsard,FoodTreatmentSatiated]"
## [46] "r_Nest[Chevroux,FoodTreatmentSatiated]"
## [47] "r_Nest[CorcellesFavres,FoodTreatmentSatiated]"
## [48] "r_Nest[Etrabloz,FoodTreatmentSatiated]"
## [49] "r_Nest[Forel,FoodTreatmentSatiated]"
## [50] "r_Nest[Franex,FoodTreatmentSatiated]"
## [51] "r_Nest[GDLV,FoodTreatmentSatiated]"
## [52] "r_Nest[Gletterens,FoodTreatmentSatiated]"
## [53] "r_Nest[Henniez,FoodTreatmentSatiated]"
## [54] "r_Nest[Jeuss,FoodTreatmentSatiated]"
## [55] "r_Nest[LesPlanches,FoodTreatmentSatiated]"
## [56] "r_Nest[Lucens,FoodTreatmentSatiated]"
## [57] "r_Nest[Lully,FoodTreatmentSatiated]"
## [58] "r_Nest[Marnand,FoodTreatmentSatiated]"
## [59] "r_Nest[Montet,FoodTreatmentSatiated]"
## [60] "r_Nest[Murist,FoodTreatmentSatiated]"
## [61] "r_Nest[Oleyes,FoodTreatmentSatiated]"
## [62] "r_Nest[Payerne,FoodTreatmentSatiated]"
## [63] "r_Nest[Rueyes,FoodTreatmentSatiated]"
## [64] "r_Nest[Seiry,FoodTreatmentSatiated]"
## [65] "r_Nest[Sevaz,FoodTreatmentSatiated]"
## [66] "r_Nest[StAubin,FoodTreatmentSatiated]"
## [67] "r_Nest[Trey,FoodTreatmentSatiated]"
## [68] "r_Nest[Yvonnand,FoodTreatmentSatiated]"
## [69] "r_Nest[AutavauxTV,SexParentMale]"
## [70] "r_Nest[Bochet,SexParentMale]"
## [71] "r_Nest[Champmartin,SexParentMale]"
## [72] "r_Nest[ChEsard,SexParentMale]"
## [73] "r_Nest[Chevroux,SexParentMale]"
## [74] "r_Nest[CorcellesFavres,SexParentMale]"
## [75] "r_Nest[Etrabloz,SexParentMale]"
## [76] "r_Nest[Forel,SexParentMale]"
## [77] "r_Nest[Franex,SexParentMale]"
## [78] "r_Nest[GDLV,SexParentMale]"
## [79] "r_Nest[Gletterens,SexParentMale]"
## [80] "r_Nest[Henniez,SexParentMale]"
## [81] "r_Nest[Jeuss,SexParentMale]"
## [82] "r_Nest[LesPlanches,SexParentMale]"
## [83] "r_Nest[Lucens,SexParentMale]"
## [84] "r_Nest[Lully,SexParentMale]"
## [85] "r_Nest[Marnand,SexParentMale]"
## [86] "r_Nest[Montet,SexParentMale]"
## [87] "r_Nest[Murist,SexParentMale]"
## [88] "r_Nest[Oleyes,SexParentMale]"
## [89] "r_Nest[Payerne,SexParentMale]"
## [90] "r_Nest[Rueyes,SexParentMale]"
## [91] "r_Nest[Seiry,SexParentMale]"
## [92] "r_Nest[Sevaz,SexParentMale]"
## [93] "r_Nest[StAubin,SexParentMale]"
## [94] "r_Nest[Trey,SexParentMale]"
## [95] "r_Nest[Yvonnand,SexParentMale]"
## [96] "r_Nest[AutavauxTV,FoodTreatmentSatiated:SexParentMale]"
## [97] "r_Nest[Bochet,FoodTreatmentSatiated:SexParentMale]"
## [98] "r_Nest[Champmartin,FoodTreatmentSatiated:SexParentMale]"
## [99] "r_Nest[ChEsard,FoodTreatmentSatiated:SexParentMale]"
## [100] "r_Nest[Chevroux,FoodTreatmentSatiated:SexParentMale]"
## [101] "r_Nest[CorcellesFavres,FoodTreatmentSatiated:SexParentMale]"
## [102] "r_Nest[Etrabloz,FoodTreatmentSatiated:SexParentMale]"
## [103] "r_Nest[Forel,FoodTreatmentSatiated:SexParentMale]"
## [104] "r_Nest[Franex,FoodTreatmentSatiated:SexParentMale]"
## [105] "r_Nest[GDLV,FoodTreatmentSatiated:SexParentMale]"
## [106] "r_Nest[Gletterens,FoodTreatmentSatiated:SexParentMale]"
## [107] "r_Nest[Henniez,FoodTreatmentSatiated:SexParentMale]"
## [108] "r_Nest[Jeuss,FoodTreatmentSatiated:SexParentMale]"
## [109] "r_Nest[LesPlanches,FoodTreatmentSatiated:SexParentMale]"
## [110] "r_Nest[Lucens,FoodTreatmentSatiated:SexParentMale]"
## [111] "r_Nest[Lully,FoodTreatmentSatiated:SexParentMale]"
## [112] "r_Nest[Marnand,FoodTreatmentSatiated:SexParentMale]"
## [113] "r_Nest[Montet,FoodTreatmentSatiated:SexParentMale]"
## [114] "r_Nest[Murist,FoodTreatmentSatiated:SexParentMale]"
## [115] "r_Nest[Oleyes,FoodTreatmentSatiated:SexParentMale]"
## [116] "r_Nest[Payerne,FoodTreatmentSatiated:SexParentMale]"
## [117] "r_Nest[Rueyes,FoodTreatmentSatiated:SexParentMale]"
## [118] "r_Nest[Seiry,FoodTreatmentSatiated:SexParentMale]"
## [119] "r_Nest[Sevaz,FoodTreatmentSatiated:SexParentMale]"
## [120] "r_Nest[StAubin,FoodTreatmentSatiated:SexParentMale]"
## [121] "r_Nest[Trey,FoodTreatmentSatiated:SexParentMale]"
## [122] "r_Nest[Yvonnand,FoodTreatmentSatiated:SexParentMale]"
## [123] "prior_Intercept"
## [124] "prior_b"
## [125] "prior_sd_Nest"
## [126] "prior_cor_Nest"
## [127] "lp__"
## [128] "accept_stat__"
## [129] "stepsize__"
## [130] "treedepth__"
## [131] "n_leapfrog__"
## [132] "divergent__"
## [133] "energy__"
owls.brm3 %>%
hypothesis("FoodTreatmentSatiated=0") %>%
plot()
owls.brm3 %>%
posterior_samples() %>%
dplyr::select(-`lp__`) %>%
pivot_longer(everything(), names_to = "key") %>%
filter(!str_detect(key, "^r")) %>%
mutate(
Type = ifelse(str_detect(key, "prior"), "Prior", "Posterior"),
## Class = ifelse(str_detect(key, 'Intercept'), 'Intercept',
## ifelse(str_detect(key, 'b'), 'b', 'sigma')),
Class = case_when(
str_detect(key, "(^b|^prior).*Intercept$") ~ "Intercept",
str_detect(key, "b_FoodTreatment.*|b_SexParent.*|prior_b") ~ "TREATMENT",
str_detect(key, "sd") ~ "sd",
str_detect(key, "^cor|prior_cor") ~ "cor",
str_detect(key, "sigma") ~ "sigma"
),
Par = str_replace(key, "b_", "")
) %>%
ggplot(aes(x = Type, y = value, color = Par)) +
stat_pointinterval(position = position_dodge(), show.legend = FALSE) +
facet_wrap(~Class, scales = "free")
In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.
There are a wide variety of tests that range from the big picture, overall chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.
The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.
available_mcmc()
## bayesplot MCMC module:
## mcmc_acf
## mcmc_acf_bar
## mcmc_areas
## mcmc_areas_data
## mcmc_areas_ridges
## mcmc_areas_ridges_data
## mcmc_combo
## mcmc_dens
## mcmc_dens_chains
## mcmc_dens_chains_data
## mcmc_dens_overlay
## mcmc_hex
## mcmc_hist
## mcmc_hist_by_chain
## mcmc_intervals
## mcmc_intervals_data
## mcmc_neff
## mcmc_neff_data
## mcmc_neff_hist
## mcmc_nuts_acceptance
## mcmc_nuts_divergence
## mcmc_nuts_energy
## mcmc_nuts_stepsize
## mcmc_nuts_treedepth
## mcmc_pairs
## mcmc_parcoord
## mcmc_parcoord_data
## mcmc_rank_hist
## mcmc_rank_overlay
## mcmc_recover_hist
## mcmc_recover_intervals
## mcmc_recover_scatter
## mcmc_rhat
## mcmc_rhat_data
## mcmc_rhat_hist
## mcmc_scatter
## mcmc_trace
## mcmc_trace_data
## mcmc_trace_highlight
## mcmc_violin
Of these, we will focus on:
pars <- owls.brm3 %>% get_variables()
pars <- pars %>%
str_extract("^b.Intercept|^b_FoodTreatment.*|^b_SexParent.*|[sS]igma|^sd.*") %>%
na.omit()
pars
## [1] "b_Intercept"
## [2] "b_FoodTreatmentSatiated"
## [3] "b_SexParentMale"
## [4] "b_FoodTreatmentSatiated:SexParentMale"
## [5] "sd_Nest__Intercept"
## [6] "sd_Nest__FoodTreatmentSatiated"
## [7] "sd_Nest__SexParentMale"
## [8] "sd_Nest__FoodTreatmentSatiated:SexParentMale"
## attr(,"na.action")
## [1] 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## [19] 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
## [37] 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
## [55] 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## [73] 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
## [91] 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
## [109] 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133
## attr(,"class")
## [1] "omit"
owls.brm3 %>% mcmc_plot(type = "trace", pars = pars)
The chains appear well mixed and very similar
owls.brm3 %>% mcmc_plot(type = "acf_bar", pars = pars)
There is no evidence of auto-correlation in the MCMC samples
owls.brm3 %>% mcmc_plot(type = "rhat_hist")
All Rhat values are below 1.05, suggesting the chains have converged.
neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
owls.brm3 %>% mcmc_plot(type = "neff_hist")
Ratios all very high.
owls.brm3 %>% mcmc_plot(type = "combo", pars = pars)
owls.brm3 %>% mcmc_plot(type = "violin", pars = pars)
The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.
Of these, we will focus on:
owls.brm3 %>% get_variables()
## [1] "b_Intercept"
## [2] "b_FoodTreatmentSatiated"
## [3] "b_SexParentMale"
## [4] "b_FoodTreatmentSatiated:SexParentMale"
## [5] "sd_Nest__Intercept"
## [6] "sd_Nest__FoodTreatmentSatiated"
## [7] "sd_Nest__SexParentMale"
## [8] "sd_Nest__FoodTreatmentSatiated:SexParentMale"
## [9] "cor_Nest__Intercept__FoodTreatmentSatiated"
## [10] "cor_Nest__Intercept__SexParentMale"
## [11] "cor_Nest__FoodTreatmentSatiated__SexParentMale"
## [12] "cor_Nest__Intercept__FoodTreatmentSatiated:SexParentMale"
## [13] "cor_Nest__FoodTreatmentSatiated__FoodTreatmentSatiated:SexParentMale"
## [14] "cor_Nest__SexParentMale__FoodTreatmentSatiated:SexParentMale"
## [15] "r_Nest[AutavauxTV,Intercept]"
## [16] "r_Nest[Bochet,Intercept]"
## [17] "r_Nest[Champmartin,Intercept]"
## [18] "r_Nest[ChEsard,Intercept]"
## [19] "r_Nest[Chevroux,Intercept]"
## [20] "r_Nest[CorcellesFavres,Intercept]"
## [21] "r_Nest[Etrabloz,Intercept]"
## [22] "r_Nest[Forel,Intercept]"
## [23] "r_Nest[Franex,Intercept]"
## [24] "r_Nest[GDLV,Intercept]"
## [25] "r_Nest[Gletterens,Intercept]"
## [26] "r_Nest[Henniez,Intercept]"
## [27] "r_Nest[Jeuss,Intercept]"
## [28] "r_Nest[LesPlanches,Intercept]"
## [29] "r_Nest[Lucens,Intercept]"
## [30] "r_Nest[Lully,Intercept]"
## [31] "r_Nest[Marnand,Intercept]"
## [32] "r_Nest[Montet,Intercept]"
## [33] "r_Nest[Murist,Intercept]"
## [34] "r_Nest[Oleyes,Intercept]"
## [35] "r_Nest[Payerne,Intercept]"
## [36] "r_Nest[Rueyes,Intercept]"
## [37] "r_Nest[Seiry,Intercept]"
## [38] "r_Nest[Sevaz,Intercept]"
## [39] "r_Nest[StAubin,Intercept]"
## [40] "r_Nest[Trey,Intercept]"
## [41] "r_Nest[Yvonnand,Intercept]"
## [42] "r_Nest[AutavauxTV,FoodTreatmentSatiated]"
## [43] "r_Nest[Bochet,FoodTreatmentSatiated]"
## [44] "r_Nest[Champmartin,FoodTreatmentSatiated]"
## [45] "r_Nest[ChEsard,FoodTreatmentSatiated]"
## [46] "r_Nest[Chevroux,FoodTreatmentSatiated]"
## [47] "r_Nest[CorcellesFavres,FoodTreatmentSatiated]"
## [48] "r_Nest[Etrabloz,FoodTreatmentSatiated]"
## [49] "r_Nest[Forel,FoodTreatmentSatiated]"
## [50] "r_Nest[Franex,FoodTreatmentSatiated]"
## [51] "r_Nest[GDLV,FoodTreatmentSatiated]"
## [52] "r_Nest[Gletterens,FoodTreatmentSatiated]"
## [53] "r_Nest[Henniez,FoodTreatmentSatiated]"
## [54] "r_Nest[Jeuss,FoodTreatmentSatiated]"
## [55] "r_Nest[LesPlanches,FoodTreatmentSatiated]"
## [56] "r_Nest[Lucens,FoodTreatmentSatiated]"
## [57] "r_Nest[Lully,FoodTreatmentSatiated]"
## [58] "r_Nest[Marnand,FoodTreatmentSatiated]"
## [59] "r_Nest[Montet,FoodTreatmentSatiated]"
## [60] "r_Nest[Murist,FoodTreatmentSatiated]"
## [61] "r_Nest[Oleyes,FoodTreatmentSatiated]"
## [62] "r_Nest[Payerne,FoodTreatmentSatiated]"
## [63] "r_Nest[Rueyes,FoodTreatmentSatiated]"
## [64] "r_Nest[Seiry,FoodTreatmentSatiated]"
## [65] "r_Nest[Sevaz,FoodTreatmentSatiated]"
## [66] "r_Nest[StAubin,FoodTreatmentSatiated]"
## [67] "r_Nest[Trey,FoodTreatmentSatiated]"
## [68] "r_Nest[Yvonnand,FoodTreatmentSatiated]"
## [69] "r_Nest[AutavauxTV,SexParentMale]"
## [70] "r_Nest[Bochet,SexParentMale]"
## [71] "r_Nest[Champmartin,SexParentMale]"
## [72] "r_Nest[ChEsard,SexParentMale]"
## [73] "r_Nest[Chevroux,SexParentMale]"
## [74] "r_Nest[CorcellesFavres,SexParentMale]"
## [75] "r_Nest[Etrabloz,SexParentMale]"
## [76] "r_Nest[Forel,SexParentMale]"
## [77] "r_Nest[Franex,SexParentMale]"
## [78] "r_Nest[GDLV,SexParentMale]"
## [79] "r_Nest[Gletterens,SexParentMale]"
## [80] "r_Nest[Henniez,SexParentMale]"
## [81] "r_Nest[Jeuss,SexParentMale]"
## [82] "r_Nest[LesPlanches,SexParentMale]"
## [83] "r_Nest[Lucens,SexParentMale]"
## [84] "r_Nest[Lully,SexParentMale]"
## [85] "r_Nest[Marnand,SexParentMale]"
## [86] "r_Nest[Montet,SexParentMale]"
## [87] "r_Nest[Murist,SexParentMale]"
## [88] "r_Nest[Oleyes,SexParentMale]"
## [89] "r_Nest[Payerne,SexParentMale]"
## [90] "r_Nest[Rueyes,SexParentMale]"
## [91] "r_Nest[Seiry,SexParentMale]"
## [92] "r_Nest[Sevaz,SexParentMale]"
## [93] "r_Nest[StAubin,SexParentMale]"
## [94] "r_Nest[Trey,SexParentMale]"
## [95] "r_Nest[Yvonnand,SexParentMale]"
## [96] "r_Nest[AutavauxTV,FoodTreatmentSatiated:SexParentMale]"
## [97] "r_Nest[Bochet,FoodTreatmentSatiated:SexParentMale]"
## [98] "r_Nest[Champmartin,FoodTreatmentSatiated:SexParentMale]"
## [99] "r_Nest[ChEsard,FoodTreatmentSatiated:SexParentMale]"
## [100] "r_Nest[Chevroux,FoodTreatmentSatiated:SexParentMale]"
## [101] "r_Nest[CorcellesFavres,FoodTreatmentSatiated:SexParentMale]"
## [102] "r_Nest[Etrabloz,FoodTreatmentSatiated:SexParentMale]"
## [103] "r_Nest[Forel,FoodTreatmentSatiated:SexParentMale]"
## [104] "r_Nest[Franex,FoodTreatmentSatiated:SexParentMale]"
## [105] "r_Nest[GDLV,FoodTreatmentSatiated:SexParentMale]"
## [106] "r_Nest[Gletterens,FoodTreatmentSatiated:SexParentMale]"
## [107] "r_Nest[Henniez,FoodTreatmentSatiated:SexParentMale]"
## [108] "r_Nest[Jeuss,FoodTreatmentSatiated:SexParentMale]"
## [109] "r_Nest[LesPlanches,FoodTreatmentSatiated:SexParentMale]"
## [110] "r_Nest[Lucens,FoodTreatmentSatiated:SexParentMale]"
## [111] "r_Nest[Lully,FoodTreatmentSatiated:SexParentMale]"
## [112] "r_Nest[Marnand,FoodTreatmentSatiated:SexParentMale]"
## [113] "r_Nest[Montet,FoodTreatmentSatiated:SexParentMale]"
## [114] "r_Nest[Murist,FoodTreatmentSatiated:SexParentMale]"
## [115] "r_Nest[Oleyes,FoodTreatmentSatiated:SexParentMale]"
## [116] "r_Nest[Payerne,FoodTreatmentSatiated:SexParentMale]"
## [117] "r_Nest[Rueyes,FoodTreatmentSatiated:SexParentMale]"
## [118] "r_Nest[Seiry,FoodTreatmentSatiated:SexParentMale]"
## [119] "r_Nest[Sevaz,FoodTreatmentSatiated:SexParentMale]"
## [120] "r_Nest[StAubin,FoodTreatmentSatiated:SexParentMale]"
## [121] "r_Nest[Trey,FoodTreatmentSatiated:SexParentMale]"
## [122] "r_Nest[Yvonnand,FoodTreatmentSatiated:SexParentMale]"
## [123] "prior_Intercept"
## [124] "prior_b"
## [125] "prior_sd_Nest"
## [126] "prior_cor_Nest"
## [127] "lp__"
## [128] "accept_stat__"
## [129] "stepsize__"
## [130] "treedepth__"
## [131] "n_leapfrog__"
## [132] "divergent__"
## [133] "energy__"
pars <- owls.brm3 %>% get_variables()
pars <- str_extract(pars, "^b_.*|^sigma$|^sd.*") %>% na.omit()
owls.brm3$fit %>%
stan_trace(pars = pars)
The chains appear well mixed and very similar
owls.brm3$fit %>%
stan_ac(pars = pars)
There is no evidence of auto-correlation in the MCMC samples
owls.brm3$fit %>% stan_rhat()
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
owls.brm3$fit %>% stan_ess()
Ratios all very high.
owls.brm3$fit %>%
stan_dens(separate_chains = TRUE, pars = pars)
The ggmean package also as a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.
Of these, we will focus on:
## owls.ggs <- owls.brm3 %>% ggs(burnin = FALSE, inc_warmup = FALSE)
## owls.ggs %>% ggs_traceplot()
The chains appear well mixed and very similar
## ggs_autocorrelation(owls.ggs)
There is no evidence of auto-correlation in the MCMC samples
## ggs_Rhat(owls.ggs)
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
## ggs_effective(owls.ggs)
Ratios all very high.
## ggs_crosscorrelation(owls.ggs)
## ggs_grb(owls.ggs)
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
available_ppc()
## bayesplot PPC module:
## ppc_bars
## ppc_bars_grouped
## ppc_boxplot
## ppc_data
## ppc_dens
## ppc_dens_overlay
## ppc_dens_overlay_grouped
## ppc_ecdf_overlay
## ppc_ecdf_overlay_grouped
## ppc_error_binned
## ppc_error_hist
## ppc_error_hist_grouped
## ppc_error_scatter
## ppc_error_scatter_avg
## ppc_error_scatter_avg_vs_x
## ppc_freqpoly
## ppc_freqpoly_grouped
## ppc_hist
## ppc_intervals
## ppc_intervals_data
## ppc_intervals_grouped
## ppc_km_overlay
## ppc_loo_intervals
## ppc_loo_pit
## ppc_loo_pit_data
## ppc_loo_pit_overlay
## ppc_loo_pit_qq
## ppc_loo_ribbon
## ppc_ribbon
## ppc_ribbon_data
## ppc_ribbon_grouped
## ppc_rootogram
## ppc_scatter
## ppc_scatter_avg
## ppc_scatter_avg_grouped
## ppc_stat
## ppc_stat_2d
## ppc_stat_freqpoly_grouped
## ppc_stat_grouped
## ppc_violin_grouped
owls.brm3 %>% pp_check(type = "dens_overlay", nsamples = 100)
The model draws appear to differ from the observed data.
owls.brm3 %>% pp_check(type = "error_scatter_avg")
This is not really interpretable
owls.brm3 %>% pp_check(group = "Nest", type = "intervals")
The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.
# library(shinystan)
# launch_shinystan(owls.brm2)
DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components yourself, we can still obtain the simulated residuals from the fitted stan model.
We need to supply:
preds <- owls.brm3 %>% posterior_predict(nsamples = 250, summary = FALSE)
owls.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = owls$NCalls,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE
)
plot(owls.resids, quantreg = TRUE)
Conclusions:
Perhaps we should specifically explore zero-inflation.
owls.resids %>% testZeroInflation()
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 3.3577, p-value < 2.2e-16
## alternative hypothesis: two.sided
Conclusions:
The data were collected at various times throughout the night. It is possible that this could lead to patterns of dependency that are not already accounted for. For example, perhaps observations that are collected at similar time of the night (within a given nest) have more similar residuals than those at very different time of the night. We can explore whether there are any temporal auto-correlation patterns.
owls.resids %>% testTemporalAutocorrelation(time = owls$ArrivalTime)
## Error in testTemporalAutocorrelation(., time = owls$ArrivalTime): testing for temporal autocorrelation requires unique time values - if you have several observations per time value, either use the recalculateResiduals function to aggregate residuals per time step, or extract the residuals from the fitted object, and plot / test each of them independently for temporally repeated subgroups (typical choices would be location / subject etc.). Note that the latter must be done by hand, outside testSpatialAutocorrelation.
owls.resid1 <- owls.resids %>% recalculateResiduals(group = interaction(owls$ArrivalTime, owls$Nest), aggregateBy = mean)
owls.resid1 %>% testTemporalAutocorrelation(time = unique(owls$ArrivalTime))
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 1.9407, p-value = 0.5966
## alternative hypothesis: true autocorrelation is not 0
Conclusions:
Conclusions:
glmer(), so we will proceed with glmmTMB() only.Zero-inflated vs hurdle models
priors <- prior(normal(1.8, 5), class = "Intercept") +
prior(normal(0, 5), class = "b") +
prior(cauchy(0, 1), class = "sd") +
prior(logistic(0, 1), class = "Intercept", dpar = "zi")
owls.form <- bf(NCalls ~ FoodTreatment * SexParent +
offset(log(BroodSize)) +
(FoodTreatment * SexParent | Nest),
zi ~ 1,
family = zero_inflated_poisson(link = "log")
)
owls.brm4 <- brm(owls.form,
data = owls,
prior = priors,
sample_prior = "yes",
iter = 5000,
warmup = 2500,
thin = 10,
chains = 3,
refresh = 0,
cores = 3
)
Model validation
preds <- owls.brm4 %>% posterior_predict(nsamples = 250, summary = FALSE)
owls.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = owls$NCalls,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE
)
plot(owls.resids, quantreg = TRUE)
owls.resids %>% testZeroInflation()
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0295, p-value = 0.792
## alternative hypothesis: two.sided
owls.resids %>% testDispersion()
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.1987, p-value = 0.048
## alternative hypothesis: two.sided
owls.resids %>% testUniformity()
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.12482, p-value = 1.569e-08
## alternative hypothesis: two-sided
owls.resids %>% testQuantiles()
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 3.413e-05
## alternative hypothesis: both
owls.resids %>% testResiduals()
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.12482, p-value = 1.569e-08
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.1987, p-value = 0.048
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 21, observations = 599, p-value < 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.009223706
## sample estimates:
## outlier frequency (expected: 0.00400667779632721 )
## 0.03505843
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.12482, p-value = 1.569e-08
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.1987, p-value = 0.048
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 21, observations = 599, p-value < 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.009223706
## sample estimates:
## outlier frequency (expected: 0.00400667779632721 )
## 0.03505843
priors <- prior(normal(1.8, 5), class = "Intercept") +
prior(normal(0, 5), class = "b") +
prior(cauchy(0, 1), class = "sd") +
prior(logistic(0, 1), class = "Intercept", dpar = "zi") +
prior(normal(0, 1), class = "b", dpar = "zi")
owls.form <- bf(NCalls ~ FoodTreatment * SexParent +
offset(log(BroodSize)) +
(FoodTreatment * SexParent | Nest),
zi ~ FoodTreatment * SexParent,
family = zero_inflated_poisson(link = "log")
)
owls.brm5 <- brm(owls.form,
data = owls,
prior = priors,
sample_prior = "yes",
iter = 5000,
warmup = 2500,
thin = 10,
chains = 3,
refresh = 0,
cores = 3
)
Model validation
preds <- owls.brm5 %>% posterior_predict(nsamples = 250, summary = FALSE)
owls.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = owls$NCalls,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE
)
plot(owls.resids, quantreg = TRUE)
owls.resids %>% testZeroInflation()
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0223, p-value = 0.832
## alternative hypothesis: two.sided
owls.resids %>% testDispersion()
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.2102, p-value = 0.04
## alternative hypothesis: two.sided
owls.resids %>% testUniformity()
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.090609, p-value = 0.000107
## alternative hypothesis: two-sided
owls.resids %>% testQuantiles()
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 8.086e-05
## alternative hypothesis: both
owls.resids %>% testResiduals()
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.090609, p-value = 0.000107
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.2102, p-value = 0.04
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 19, observations = 599, p-value < 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01168614
## sample estimates:
## outlier frequency (expected: 0.00452420701168614 )
## 0.03171953
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.090609, p-value = 0.000107
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 1.2102, p-value = 0.04
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 19, observations = 599, p-value < 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01168614
## sample estimates:
## outlier frequency (expected: 0.00452420701168614 )
## 0.03171953
priors <- prior(normal(1.8, 5), class = "Intercept") +
prior(normal(0, 5), class = "b") +
prior(cauchy(0, 1), class = "sd") +
prior(logistic(0, 1), class = "Intercept", dpar = "zi") +
prior(gamma(0.01, 0.01), class = "shape")
owls.form <- bf(NCalls ~ FoodTreatment * SexParent +
offset(log(BroodSize)) +
(FoodTreatment * SexParent | Nest),
zi ~ 1,
family = zero_inflated_negbinomial(link = "log")
)
owls.brm6 <- brm(owls.form,
data = owls,
prior = priors,
sample_prior = "yes",
iter = 5000,
warmup = 2500,
thin = 10,
chains = 3,
refresh = 0,
cores = 3
)
Model validation
preds <- owls.brm6 %>% posterior_predict(nsamples = 250, summary = FALSE)
owls.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = owls$NCalls,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE
)
plot(owls.resids, quantreg = TRUE)
owls.resids %>% testZeroInflation()
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.048, p-value = 0.592
## alternative hypothesis: two.sided
owls.resids %>% testDispersion()
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.83884, p-value = 0.032
## alternative hypothesis: two.sided
owls.resids %>% testUniformity()
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.033463, p-value = 0.5136
## alternative hypothesis: two-sided
owls.resids %>% testQuantiles()
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.00197
## alternative hypothesis: both
owls.resids %>% testResiduals()
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.033463, p-value = 0.5136
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.83884, p-value = 0.032
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 599, p-value = 0.18
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01335559
## sample estimates:
## outlier frequency (expected: 0.00534223706176962 )
## 0
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.033463, p-value = 0.5136
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.83884, p-value = 0.032
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 599, p-value = 0.18
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01335559
## sample estimates:
## outlier frequency (expected: 0.00534223706176962 )
## 0
priors <- prior(normal(1.8, 5), class = "Intercept") +
prior(normal(0, 5), class = "b") +
prior(cauchy(0, 1), class = "sd") +
prior(logistic(0, 1), class = "Intercept", dpar = "zi") +
prior(normal(0, 1), class = "b", dpar = "zi") +
prior(gamma(0.01, 0.01), class = "shape")
owls.form <- bf(NCalls ~ FoodTreatment * SexParent +
offset(log(BroodSize)) +
(FoodTreatment * SexParent | Nest),
zi ~ FoodTreatment * SexParent,
family = zero_inflated_negbinomial(link = "log")
)
owls.brm7 <- brm(owls.form,
data = owls,
prior = priors,
sample_prior = "yes",
iter = 5000,
warmup = 2500,
thin = 10,
chains = 3,
refresh = 0,
cores = 3
)
Model validation
preds <- owls.brm7 %>% posterior_predict(nsamples = 250, summary = FALSE)
owls.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = owls$NCalls,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = TRUE
)
plot(owls.resids, quantreg = TRUE)
owls.resids %>% testZeroInflation()
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0258, p-value = 0.824
## alternative hypothesis: two.sided
owls.resids %>% testDispersion()
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.84071, p-value = 0.016
## alternative hypothesis: two.sided
owls.resids %>% testUniformity()
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.044661, p-value = 0.1832
## alternative hypothesis: two-sided
owls.resids %>% testQuantiles()
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.006258
## alternative hypothesis: both
owls.resids %>% testResiduals()
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.044661, p-value = 0.1832
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.84071, p-value = 0.016
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 599, p-value = 0.08
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01335559
## sample estimates:
## outlier frequency (expected: 0.00529215358931553 )
## 0
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.044661, p-value = 0.1832
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.84071, p-value = 0.016
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 599, p-value = 0.08
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01335559
## sample estimates:
## outlier frequency (expected: 0.00529215358931553 )
## 0
(l.1 <- loo(owls.brm3))
##
## Computed from 750 by 599 log-likelihood matrix
##
## Estimate SE
## elpd_loo -2430.8 75.5
## p_loo 301.3 20.8
## looic 4861.6 150.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 507 84.6% 81
## (0.5, 0.7] (ok) 49 8.2% 51
## (0.7, 1] (bad) 37 6.2% 19
## (1, Inf) (very bad) 6 1.0% 1
## See help('pareto-k-diagnostic') for details.
(l.2 <- loo(owls.brm4))
##
## Computed from 750 by 599 log-likelihood matrix
##
## Estimate SE
## elpd_loo -1930.9 53.3
## p_loo 185.3 13.4
## looic 3861.8 106.6
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 528 88.1% 88
## (0.5, 0.7] (ok) 42 7.0% 43
## (0.7, 1] (bad) 25 4.2% 19
## (1, Inf) (very bad) 4 0.7% 6
## See help('pareto-k-diagnostic') for details.
(l.3 <- loo(owls.brm5))
##
## Computed from 750 by 599 log-likelihood matrix
##
## Estimate SE
## elpd_loo -1920.0 53.8
## p_loo 184.4 12.9
## looic 3840.0 107.7
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 526 87.8% 62
## (0.5, 0.7] (ok) 48 8.0% 41
## (0.7, 1] (bad) 18 3.0% 15
## (1, Inf) (very bad) 7 1.2% 4
## See help('pareto-k-diagnostic') for details.
(l.4 <- loo(owls.brm6))
##
## Computed from 750 by 599 log-likelihood matrix
##
## Estimate SE
## elpd_loo -1665.6 28.4
## p_loo 50.7 4.8
## looic 3331.1 56.8
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 582 97.2% 150
## (0.5, 0.7] (ok) 13 2.2% 82
## (0.7, 1] (bad) 3 0.5% 23
## (1, Inf) (very bad) 1 0.2% 12
## See help('pareto-k-diagnostic') for details.
(l.5 <- loo(owls.brm7))
##
## Computed from 750 by 599 log-likelihood matrix
##
## Estimate SE
## elpd_loo -1656.1 29.0
## p_loo 50.7 4.1
## looic 3312.3 57.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 575 96.0% 167
## (0.5, 0.7] (ok) 17 2.8% 75
## (0.7, 1] (bad) 6 1.0% 64
## (1, Inf) (very bad) 1 0.2% 9
## See help('pareto-k-diagnostic') for details.
loo_compare(l.1, l.2, l.3, l.4, l.5)
## elpd_diff se_diff
## owls.brm7 0.0 0.0
## owls.brm6 -9.4 5.0
## owls.brm5 -263.8 34.7
## owls.brm4 -274.8 35.0
## owls.brm3 -774.6 66.6
owls.brm7 %>%
conditional_effects("FoodTreatment:SexParent") %>%
plot(points = TRUE, jitter = c(0.25, 0))
These predictions appear to be based on the mean BroodSize of approximately 4.39. That is, the predictions are for a nest, not per chick. There does not appear to be a way to indicate the offset value.
owls.brm7 %>%
ggpredict(~ FoodTreatment * SexParent) %>%
plot(add.data = TRUE, jitter = c(0.25, 0))
These predictions appear to be based on the mean BroodSize of approximately 4.39. That is, the predictions are for a nest, not per chick. There does not appear to be a way to indicate the offset value.
ggemmeans() can accommodate the offset correctly. There are two sensible choices:
off <- owls %>% summarize(Mean = mean(BroodSize))
as.numeric(off)
## [1] 4.392321
owls.brm7 %>%
ggemmeans(~ FoodTreatment * SexParent, offset = off$Mean) %>%
plot()
owls.brm7 %>%
ggemmeans(~ FoodTreatment * SexParent, offset = log(off$Mean)) %>%
plot()
owls.brm7 %>%
ggemmeans(~ FoodTreatment * SexParent, offset = log(1)) %>%
plot()
brms captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.
The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).
owls.brm7 %>% summary()
## Family: zero_inflated_negbinomial
## Links: mu = log; shape = identity; zi = logit
## Formula: NCalls ~ FoodTreatment * SexParent + offset(log(BroodSize)) + (FoodTreatment * SexParent | Nest)
## zi ~ FoodTreatment * SexParent
## Data: owls (Number of observations: 599)
## Samples: 3 chains, each with iter = 5000; warmup = 2500; thin = 10;
## total post-warmup samples = 750
##
## Group-Level Effects:
## ~Nest (Number of levels: 27)
## Estimate
## sd(Intercept) 0.30
## sd(FoodTreatmentSatiated) 1.17
## sd(SexParentMale) 0.17
## sd(FoodTreatmentSatiated:SexParentMale) 0.36
## cor(Intercept,FoodTreatmentSatiated) 0.09
## cor(Intercept,SexParentMale) -0.17
## cor(FoodTreatmentSatiated,SexParentMale) -0.23
## cor(Intercept,FoodTreatmentSatiated:SexParentMale) -0.10
## cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale) -0.09
## cor(SexParentMale,FoodTreatmentSatiated:SexParentMale) -0.03
## Est.Error
## sd(Intercept) 0.09
## sd(FoodTreatmentSatiated) 0.29
## sd(SexParentMale) 0.12
## sd(FoodTreatmentSatiated:SexParentMale) 0.23
## cor(Intercept,FoodTreatmentSatiated) 0.29
## cor(Intercept,SexParentMale) 0.40
## cor(FoodTreatmentSatiated,SexParentMale) 0.41
## cor(Intercept,FoodTreatmentSatiated:SexParentMale) 0.40
## cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale) 0.43
## cor(SexParentMale,FoodTreatmentSatiated:SexParentMale) 0.43
## l-95% CI
## sd(Intercept) 0.14
## sd(FoodTreatmentSatiated) 0.66
## sd(SexParentMale) 0.01
## sd(FoodTreatmentSatiated:SexParentMale) 0.02
## cor(Intercept,FoodTreatmentSatiated) -0.45
## cor(Intercept,SexParentMale) -0.82
## cor(FoodTreatmentSatiated,SexParentMale) -0.86
## cor(Intercept,FoodTreatmentSatiated:SexParentMale) -0.82
## cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale) -0.80
## cor(SexParentMale,FoodTreatmentSatiated:SexParentMale) -0.81
## u-95% CI Rhat
## sd(Intercept) 0.50 1.00
## sd(FoodTreatmentSatiated) 1.77 1.00
## sd(SexParentMale) 0.43 1.00
## sd(FoodTreatmentSatiated:SexParentMale) 0.87 1.00
## cor(Intercept,FoodTreatmentSatiated) 0.62 1.00
## cor(Intercept,SexParentMale) 0.71 1.00
## cor(FoodTreatmentSatiated,SexParentMale) 0.69 1.00
## cor(Intercept,FoodTreatmentSatiated:SexParentMale) 0.63 1.00
## cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale) 0.75 1.00
## cor(SexParentMale,FoodTreatmentSatiated:SexParentMale) 0.74 1.00
## Bulk_ESS
## sd(Intercept) 762
## sd(FoodTreatmentSatiated) 750
## sd(SexParentMale) 586
## sd(FoodTreatmentSatiated:SexParentMale) 799
## cor(Intercept,FoodTreatmentSatiated) 661
## cor(Intercept,SexParentMale) 754
## cor(FoodTreatmentSatiated,SexParentMale) 853
## cor(Intercept,FoodTreatmentSatiated:SexParentMale) 598
## cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale) 929
## cor(SexParentMale,FoodTreatmentSatiated:SexParentMale) 547
## Tail_ESS
## sd(Intercept) 759
## sd(FoodTreatmentSatiated) 753
## sd(SexParentMale) 754
## sd(FoodTreatmentSatiated:SexParentMale) 655
## cor(Intercept,FoodTreatmentSatiated) 695
## cor(Intercept,SexParentMale) 628
## cor(FoodTreatmentSatiated,SexParentMale) 797
## cor(Intercept,FoodTreatmentSatiated:SexParentMale) 670
## cor(FoodTreatmentSatiated,FoodTreatmentSatiated:SexParentMale) 682
## cor(SexParentMale,FoodTreatmentSatiated:SexParentMale) 512
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## Intercept 0.84 0.10 0.64 1.06
## zi_Intercept -1.61 0.25 -2.12 -1.15
## FoodTreatmentSatiated -0.77 0.29 -1.37 -0.21
## SexParentMale -0.10 0.11 -0.33 0.10
## FoodTreatmentSatiated:SexParentMale 0.13 0.20 -0.26 0.52
## zi_FoodTreatmentSatiated 0.81 0.34 0.14 1.48
## zi_SexParentMale -0.78 0.34 -1.53 -0.12
## zi_FoodTreatmentSatiated:SexParentMale 0.60 0.44 -0.26 1.52
## Rhat Bulk_ESS Tail_ESS
## Intercept 1.00 810 759
## zi_Intercept 1.00 676 691
## FoodTreatmentSatiated 1.00 634 631
## SexParentMale 1.00 781 682
## FoodTreatmentSatiated:SexParentMale 1.01 662 649
## zi_FoodTreatmentSatiated 1.00 655 630
## zi_SexParentMale 1.00 793 756
## zi_FoodTreatmentSatiated:SexParentMale 1.00 757 799
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 2.75 0.32 2.13 3.34 1.01 812 717
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Conclusions:
owls.brm7$fit %>%
tidyMCMC(
estimate.method = "median",
conf.int = TRUE, conf.method = "HPDinterval",
rhat = TRUE, ess = TRUE
)
Conclusions:
see above
Due to the presence of a log transform in the predictor, it is better to use the regex version.
owls.brm7 %>% get_variables()
## [1] "b_Intercept"
## [2] "b_zi_Intercept"
## [3] "b_FoodTreatmentSatiated"
## [4] "b_SexParentMale"
## [5] "b_FoodTreatmentSatiated:SexParentMale"
## [6] "b_zi_FoodTreatmentSatiated"
## [7] "b_zi_SexParentMale"
## [8] "b_zi_FoodTreatmentSatiated:SexParentMale"
## [9] "sd_Nest__Intercept"
## [10] "sd_Nest__FoodTreatmentSatiated"
## [11] "sd_Nest__SexParentMale"
## [12] "sd_Nest__FoodTreatmentSatiated:SexParentMale"
## [13] "cor_Nest__Intercept__FoodTreatmentSatiated"
## [14] "cor_Nest__Intercept__SexParentMale"
## [15] "cor_Nest__FoodTreatmentSatiated__SexParentMale"
## [16] "cor_Nest__Intercept__FoodTreatmentSatiated:SexParentMale"
## [17] "cor_Nest__FoodTreatmentSatiated__FoodTreatmentSatiated:SexParentMale"
## [18] "cor_Nest__SexParentMale__FoodTreatmentSatiated:SexParentMale"
## [19] "shape"
## [20] "r_Nest[AutavauxTV,Intercept]"
## [21] "r_Nest[Bochet,Intercept]"
## [22] "r_Nest[Champmartin,Intercept]"
## [23] "r_Nest[ChEsard,Intercept]"
## [24] "r_Nest[Chevroux,Intercept]"
## [25] "r_Nest[CorcellesFavres,Intercept]"
## [26] "r_Nest[Etrabloz,Intercept]"
## [27] "r_Nest[Forel,Intercept]"
## [28] "r_Nest[Franex,Intercept]"
## [29] "r_Nest[GDLV,Intercept]"
## [30] "r_Nest[Gletterens,Intercept]"
## [31] "r_Nest[Henniez,Intercept]"
## [32] "r_Nest[Jeuss,Intercept]"
## [33] "r_Nest[LesPlanches,Intercept]"
## [34] "r_Nest[Lucens,Intercept]"
## [35] "r_Nest[Lully,Intercept]"
## [36] "r_Nest[Marnand,Intercept]"
## [37] "r_Nest[Montet,Intercept]"
## [38] "r_Nest[Murist,Intercept]"
## [39] "r_Nest[Oleyes,Intercept]"
## [40] "r_Nest[Payerne,Intercept]"
## [41] "r_Nest[Rueyes,Intercept]"
## [42] "r_Nest[Seiry,Intercept]"
## [43] "r_Nest[Sevaz,Intercept]"
## [44] "r_Nest[StAubin,Intercept]"
## [45] "r_Nest[Trey,Intercept]"
## [46] "r_Nest[Yvonnand,Intercept]"
## [47] "r_Nest[AutavauxTV,FoodTreatmentSatiated]"
## [48] "r_Nest[Bochet,FoodTreatmentSatiated]"
## [49] "r_Nest[Champmartin,FoodTreatmentSatiated]"
## [50] "r_Nest[ChEsard,FoodTreatmentSatiated]"
## [51] "r_Nest[Chevroux,FoodTreatmentSatiated]"
## [52] "r_Nest[CorcellesFavres,FoodTreatmentSatiated]"
## [53] "r_Nest[Etrabloz,FoodTreatmentSatiated]"
## [54] "r_Nest[Forel,FoodTreatmentSatiated]"
## [55] "r_Nest[Franex,FoodTreatmentSatiated]"
## [56] "r_Nest[GDLV,FoodTreatmentSatiated]"
## [57] "r_Nest[Gletterens,FoodTreatmentSatiated]"
## [58] "r_Nest[Henniez,FoodTreatmentSatiated]"
## [59] "r_Nest[Jeuss,FoodTreatmentSatiated]"
## [60] "r_Nest[LesPlanches,FoodTreatmentSatiated]"
## [61] "r_Nest[Lucens,FoodTreatmentSatiated]"
## [62] "r_Nest[Lully,FoodTreatmentSatiated]"
## [63] "r_Nest[Marnand,FoodTreatmentSatiated]"
## [64] "r_Nest[Montet,FoodTreatmentSatiated]"
## [65] "r_Nest[Murist,FoodTreatmentSatiated]"
## [66] "r_Nest[Oleyes,FoodTreatmentSatiated]"
## [67] "r_Nest[Payerne,FoodTreatmentSatiated]"
## [68] "r_Nest[Rueyes,FoodTreatmentSatiated]"
## [69] "r_Nest[Seiry,FoodTreatmentSatiated]"
## [70] "r_Nest[Sevaz,FoodTreatmentSatiated]"
## [71] "r_Nest[StAubin,FoodTreatmentSatiated]"
## [72] "r_Nest[Trey,FoodTreatmentSatiated]"
## [73] "r_Nest[Yvonnand,FoodTreatmentSatiated]"
## [74] "r_Nest[AutavauxTV,SexParentMale]"
## [75] "r_Nest[Bochet,SexParentMale]"
## [76] "r_Nest[Champmartin,SexParentMale]"
## [77] "r_Nest[ChEsard,SexParentMale]"
## [78] "r_Nest[Chevroux,SexParentMale]"
## [79] "r_Nest[CorcellesFavres,SexParentMale]"
## [80] "r_Nest[Etrabloz,SexParentMale]"
## [81] "r_Nest[Forel,SexParentMale]"
## [82] "r_Nest[Franex,SexParentMale]"
## [83] "r_Nest[GDLV,SexParentMale]"
## [84] "r_Nest[Gletterens,SexParentMale]"
## [85] "r_Nest[Henniez,SexParentMale]"
## [86] "r_Nest[Jeuss,SexParentMale]"
## [87] "r_Nest[LesPlanches,SexParentMale]"
## [88] "r_Nest[Lucens,SexParentMale]"
## [89] "r_Nest[Lully,SexParentMale]"
## [90] "r_Nest[Marnand,SexParentMale]"
## [91] "r_Nest[Montet,SexParentMale]"
## [92] "r_Nest[Murist,SexParentMale]"
## [93] "r_Nest[Oleyes,SexParentMale]"
## [94] "r_Nest[Payerne,SexParentMale]"
## [95] "r_Nest[Rueyes,SexParentMale]"
## [96] "r_Nest[Seiry,SexParentMale]"
## [97] "r_Nest[Sevaz,SexParentMale]"
## [98] "r_Nest[StAubin,SexParentMale]"
## [99] "r_Nest[Trey,SexParentMale]"
## [100] "r_Nest[Yvonnand,SexParentMale]"
## [101] "r_Nest[AutavauxTV,FoodTreatmentSatiated:SexParentMale]"
## [102] "r_Nest[Bochet,FoodTreatmentSatiated:SexParentMale]"
## [103] "r_Nest[Champmartin,FoodTreatmentSatiated:SexParentMale]"
## [104] "r_Nest[ChEsard,FoodTreatmentSatiated:SexParentMale]"
## [105] "r_Nest[Chevroux,FoodTreatmentSatiated:SexParentMale]"
## [106] "r_Nest[CorcellesFavres,FoodTreatmentSatiated:SexParentMale]"
## [107] "r_Nest[Etrabloz,FoodTreatmentSatiated:SexParentMale]"
## [108] "r_Nest[Forel,FoodTreatmentSatiated:SexParentMale]"
## [109] "r_Nest[Franex,FoodTreatmentSatiated:SexParentMale]"
## [110] "r_Nest[GDLV,FoodTreatmentSatiated:SexParentMale]"
## [111] "r_Nest[Gletterens,FoodTreatmentSatiated:SexParentMale]"
## [112] "r_Nest[Henniez,FoodTreatmentSatiated:SexParentMale]"
## [113] "r_Nest[Jeuss,FoodTreatmentSatiated:SexParentMale]"
## [114] "r_Nest[LesPlanches,FoodTreatmentSatiated:SexParentMale]"
## [115] "r_Nest[Lucens,FoodTreatmentSatiated:SexParentMale]"
## [116] "r_Nest[Lully,FoodTreatmentSatiated:SexParentMale]"
## [117] "r_Nest[Marnand,FoodTreatmentSatiated:SexParentMale]"
## [118] "r_Nest[Montet,FoodTreatmentSatiated:SexParentMale]"
## [119] "r_Nest[Murist,FoodTreatmentSatiated:SexParentMale]"
## [120] "r_Nest[Oleyes,FoodTreatmentSatiated:SexParentMale]"
## [121] "r_Nest[Payerne,FoodTreatmentSatiated:SexParentMale]"
## [122] "r_Nest[Rueyes,FoodTreatmentSatiated:SexParentMale]"
## [123] "r_Nest[Seiry,FoodTreatmentSatiated:SexParentMale]"
## [124] "r_Nest[Sevaz,FoodTreatmentSatiated:SexParentMale]"
## [125] "r_Nest[StAubin,FoodTreatmentSatiated:SexParentMale]"
## [126] "r_Nest[Trey,FoodTreatmentSatiated:SexParentMale]"
## [127] "r_Nest[Yvonnand,FoodTreatmentSatiated:SexParentMale]"
## [128] "prior_Intercept"
## [129] "prior_b"
## [130] "prior_shape"
## [131] "prior_b_zi"
## [132] "prior_Intercept_zi"
## [133] "prior_sd_Nest"
## [134] "prior_cor_Nest"
## [135] "lp__"
## [136] "accept_stat__"
## [137] "stepsize__"
## [138] "treedepth__"
## [139] "n_leapfrog__"
## [140] "divergent__"
## [141] "energy__"
owls.draw <- owls.brm7 %>%
gather_draws(`b.Intercept.*|b_FoodTreatment.*|b_SexParent.*`, regex = TRUE)
owls.draw
We can then summarise this
owls.draw %>% median_hdci()
owls.brm7 %>%
gather_draws(`b_Intercept.*|b_FoodTreatment.*|b_SexParent.*`, regex = TRUE) %>%
ggplot() +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_slab(aes(
x = .value, y = .variable,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()
))
), color = "black") +
scale_fill_brewer("Interval", direction = -1, na.translate = FALSE)
owls.brm7 %>%
gather_draws(`.Intercept.*|b_FoodTreatment.*|b_SexParent.*`, regex = TRUE) %>%
ggplot() +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_halfeye(aes(x = .value, y = .variable)) +
theme_classic()
owls.brm7$fit %>% plot(type = "intervals")
This is purely a graphical depiction on the posteriors.
owls.brm7 %>% tidy_draws()
owls.brm7 %>% spread_draws(`.*Intercept.*|b_FoodTreatment.*|b_SexParent.*`, regex = TRUE)
owls.brm7 %>%
posterior_samples() %>%
as_tibble()
owls.brm7 %>%
bayes_R2(re.form = NA, summary = FALSE) %>%
median_hdci()
owls.brm7 %>%
bayes_R2(re.form = ~ (1 | Nest), summary = FALSE) %>%
median_hdci()
owls.brm7 %>%
bayes_R2(re.form = ~ (FoodTreatment * SexParent | Nest), summary = FALSE) %>%
median_hdci()
owls.brm7 %>%
emmeans(~FoodTreatment, type = "response")
## FoodTreatment prob lower.HPD upper.HPD
## Deprived 2.2 1.865 2.59
## Satiated 1.1 0.516 1.71
##
## Point estimate displayed: median
## Results are back-transformed from the log scale
## HPD interval probability: 0.95
## FoodTreatment prob lower.HPD upper.HPD
## Deprived 2.2 1.865 2.59
## Satiated 1.1 0.516 1.71
##
## Point estimate displayed: median
## Results are back-transformed from the log scale
## HPD interval probability: 0.95